/*
 * Copyright (c) 2011 Robert Kirchgessner <rkirchge@gmail.com>
 *                    Carlo Pascoe <carlo.pascoe@gmail.com>
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 * 
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include "motifsearch.h"
#include "graph.h"
#include "hashgraph.h"
#include "motif.h"

#define BIAS6 63

void
motif_template_generate(char *g6str, uint32_t *size, struct graph *motif)
{
  uint32_t *adj_mat;
//  char g6str[10];
  uint32_t nV;
  uint32_t bit_count, char_count, col_count, row_count;

//  fscanf(g6fp, "%10s", g6str);
  for(char_count = 0; g6str[char_count] != '\0'; char_count++) {
    if(g6str[char_count] == '\n') {
      g6str[char_count] = '\0';
    }
  }

  //G6 string represented as N(x) R(x)
  assert(g6str[0] != 126); //Dont handle this case!

  nV = g6str[0] - BIAS6;
  if(size) {
    *size = nV;
  }
  
  adj_mat = (uint32_t *)malloc(nV*nV*sizeof(uint32_t));
  memset(adj_mat, 0, sizeof(uint32_t)*nV*nV);
  assert(adj_mat != NULL);
  
  char_count = 1; /*Set the character counter to the first character in the string*/
  bit_count = 5;
  col_count = 1;
  uint32_t adj_val;

  while(g6str[char_count] != '\0') {
    //Iterate column wise
    if(col_count == nV) {
      break;
    }

    for(row_count = 0; row_count < col_count; row_count++) {
      adj_val = g6str[char_count]-BIAS6;
      adj_mat[COORD(row_count, col_count, nV)] = ((adj_val>>bit_count)&0x1);
      adj_mat[COORD(col_count, row_count, nV)] = ((adj_val>>bit_count)&0x1);

      if(bit_count == 0) {
        bit_count = 5;
	char_count ++;
      }
      else {
        bit_count --; 
      }
    }
    col_count++;
  }

  //Build the graph from the adjacency matrix
  graph_from_adj(adj_mat, nV, motif);  

  free(adj_mat);
}
#if 0
void
motif_template_generate(FILE *g6fp, uint32_t *size, struct graph *motif)
{
  uint32_t *adj_mat;
  char g6str[10];
  uint32_t nV;
  uint32_t bit_count, char_count, col_count, row_count;

  fscanf(g6fp, "%10s", g6str);
  for(char_count = 0; g6str[char_count] != '\0'; char_count++) {
    if(g6str[char_count] == '\n') {
      g6str[char_count] = '\0';
    }
  }

  //G6 string represented as N(x) R(x)
  assert(g6str[0] != 126); //Dont handle this case!

  nV = g6str[0] - BIAS6;
  if(size) {
    *size = nV;
  }
  
  adj_mat = (uint32_t *)malloc(nV*nV*sizeof(uint32_t));
  memset(adj_mat, 0, sizeof(uint32_t)*nV*nV);
  assert(adj_mat != NULL);
  
  char_count = 1; /*Set the character counter to the first character in the string*/
  bit_count = 5;
  col_count = 1;
  uint32_t adj_val;

  while(g6str[char_count] != '\0') {
    //Iterate column wise
    if(col_count == nV) {
      break;
    }

    for(row_count = 0; row_count < col_count; row_count++) {
      adj_val = g6str[char_count]-BIAS6;
      adj_mat[COORD(row_count, col_count, nV)] = ((adj_val>>bit_count)&0x1);
      adj_mat[COORD(col_count, row_count, nV)] = ((adj_val>>bit_count)&0x1);

      if(bit_count == 0) {
        bit_count = 5;
	char_count ++;
      }
      else {
        bit_count --; 
      }
    }
    col_count++;
  }

  //Build the graph from the adjacency matrix
  graph_from_adj(adj_mat, nV, motif);  

  free(adj_mat);
}
#endif
uint32_t
motif_search_auto(struct hash_graph *hg, uint32_t threshold, uint32_t k, struct motif_pattern *mp)
{
  uint32_t *sim_mat, *census_adj, *tmp_adj;
  uint32_t x, y, xp, yp, matched;

  //Allocated similarity matrix
  sim_mat = (uint32_t *)malloc(hg->size*hg->size*sizeof(uint32_t));
  assert(sim_mat != NULL);
  memset(sim_mat,0,hg->size*hg->size*sizeof(uint32_t));

  census_adj = (uint32_t *)malloc(k*k*sizeof(uint32_t));
  assert(census_adj != NULL);
  memset(census_adj, 0, k*k*sizeof(uint32_t));


  tmp_adj = (uint32_t *)malloc(k*k*sizeof(uint32_t));
  assert(tmp_adj != NULL);
  memset(tmp_adj, 0, k*k*sizeof(uint32_t));

  //Calculate Scoring matrix
  for(x=0; x<hg->size; x++) {
    for(y=x; y<hg->size; y++) {
      sim_mat[COORD(x,y,hg->size)] = hash_dist(hg->hash_array[x], hg->hash_array[y]);
      sim_mat[COORD(y,x,hg->size)] = hash_dist(hg->hash_array[x], hg->hash_array[y]);
    }
  }

 for(x=0; x<hg->size; x++) {
    matched = 1;
    overlay_subgraph_adj(hg->subgraph_array[x], census_adj, k);
    for(y=x; y<hg->size; y++) {
      if(sim_mat[COORD(x,y,hg->size)] < threshold) {
        matched ++;
        overlay_subgraph_adj(hg->subgraph_array[x], tmp_adj, k);

	for(xp=0; xp<k; xp++) {
	  for(yp=0;yp<k; yp++) {
	    census_adj[COORD(xp,yp,k)] += tmp_adj[COORD(xp,yp,k)];
	  }
	}
      }
    }

    for(xp=0; xp<k; xp++) {
      for(yp=0;yp<k; yp++) {
        census_adj[COORD(xp,yp,k)] /= matched;
      }
    }

    if(matched > 0) {
      motif_pattern_list_add(mp, census_adj, k, matched);
    }
  }

  free(tmp_adj);
  free(census_adj);
  free(sim_mat);
  return 0;
}

void 
edge_sort(uint32_t *Sc, int *Order, uint32_t feature_size)
{

    uint32_t dtmp;
//  int otmp;
    int itmp;
    int i, j;
    for (i = 0; i < feature_size; i++)
    {
        dtmp = Sc[0];
//    otmp = Order[0];
        itmp = 0;
        for (j = 1; j < feature_size-i; j++)
        {
            if (Sc[j] < dtmp)
            {
                dtmp = Sc[j];
//        otmp = Order[j];
                itmp = j;
            }
        }
        Sc[itmp] = Sc[feature_size-i-1];
//    Order[itmp] = Order[feature_size-i-1];
        Sc[feature_size-i-1] = dtmp;
//    Order[feature_size-i-1] = otmp;
    }

}

void 
hash_func(struct subgraph *sg, hash_type *hash, uint32_t feature_size)
{
    uint32_t *adj = malloc(feature_size*feature_size*sizeof(uint32_t));
    uint32_t *num_edges = malloc(feature_size*sizeof(uint32_t));
    overlay_subgraph_adj(sg, adj,feature_size);

    int i, j;
    uint32_t th = 0;

    for ( i = 0; i < feature_size; i++)
    {
        num_edges[i] = 0;
        for ( j = 0; j < feature_size; j++)
        {
            num_edges[i] += adj[feature_size*i+j];
        }
    }

    edge_sort(num_edges, NULL, feature_size);
    for ( i = feature_size-1; i >= 0; i--)
    {
        th = (th << 4) | num_edges[i];
    }
    *hash = th;

    free(num_edges);
    free(adj);
}

void 
free_motif_list(struct motif_list_head * mlh)
{
//will only free the motif_list data. will still nead to free the motif_list_head if dynamically allocated
    struct motif_list * tmp1;
    struct motif_list * tmp2;
    if (mlh->next != NULL)
    {
        tmp1 = mlh->next;
        while (tmp1 != NULL)
        {
            tmp2 = tmp1->next;
            free(tmp1->hg_node_map);
            free(tmp1->hgq_node_map);
            free(tmp1);
            tmp1 = tmp2;
        }
    }
}

void 
motif_align(struct hash_graph *hg, struct hash_graph *hgq, struct motif_list_head *ml, uint32_t threshold, uint32_t feature_size)
{
    assert(hgq->size > 1);
    assert(hgq->size < 5);
    int x, y, z, z1, i, tmp, tmpm, mx, my, mz, mz1, list_cnt, c1, c2, c3, c4, c5, c6;
    uint32_t *hscrs = (uint32_t *)malloc(hg->size*hgq->size*sizeof(uint32_t));

    for (x=0;x<hg->size;x++)
    {
        for (y=0;y < hgq->size;y++)
        {
            hscrs[hgq->size*x + y] = 0;
            for (i = 0; i < feature_size; i++)
            {
                hscrs[hgq->size*x + y]+= abs((15 & (hg->hash_array[x] >> 4*i)) - (15 & (hgq->hash_array[y] >> 4*i)));
            }
        }
    }

    struct motif_list *head = NULL;
    struct motif_list *tmp1 = NULL;
    struct motif_list *tmp2 = NULL;
    struct motif_list *tmp3 = NULL;
    struct motif_list *tmp4 = NULL;
    struct motif_list *tmp5 = NULL;
    struct motif_list *tmp6 = NULL;

    head = (struct motif_list *)malloc(sizeof(struct motif_list));
    tmp1 = head;
    list_cnt = 0;
    if (hgq->size == 2)
    {
        for (x=0; x < hg->size; x++)
        {
            tmpm = threshold;
            for (y=0; y < hg->size; y++)
            {
               if(x != y){
                if (hg->hash_adj[hg->size*x+y] > 0)
                {
                    tmp = (abs(hg->hash_adj[hg->size*x+y] - hgq->hash_adj[1]) + 1)*(hscrs[hgq->size*x+0] + hscrs[hgq->size*y+1]);
                    if (tmp < tmpm)
                    {
                        tmpm = tmp;
                        mx = x;
                        my = y;
                    }
                }
}
            }
            if (tmpm < threshold)
            {
                list_cnt++;
                tmp1->hg_node_map = (uint32_t *)malloc(2*sizeof(uint32_t));
                tmp1->hgq_node_map = (uint32_t *)malloc(2*sizeof(uint32_t));
                tmp1->score = (double)tmpm;
                tmp1->hg_node_map[0] = mx;
                tmp1->hg_node_map[1] = my;
                tmp1->hgq_node_map[0] = 0;
                tmp1->hgq_node_map[1] = 1;
                tmp1->next = (struct motif_list *)malloc(sizeof(struct motif_list));
                tmp2 = tmp1;
                tmp1 = tmp1->next;
            }
        }
if(list_cnt == 0){
  free(head);
  head = NULL;
}else{
        free(tmp2->next);
        tmp2->next = NULL;
}
    }
    else if (hgq->size == 3)
    {
        c1 = 0;
        c2 = 0;
        c3 = 0;
        if (hgq->hash_adj[1] > 0)c1 = 1;
        if (hgq->hash_adj[2] > 0)c2 = 1;
        if (hgq->hash_adj[5] > 0)c3 = 1;
        for (x=0; x < hg->size; x++)
        {
            for (y=0; y < hg->size; y++)
            {
                if (x != y)
                {
                    tmpm = threshold;
                    for (z=0; z < hg->size; z++)
                    {
                        if (x != z && y != z)
                        {
                            tmp = 0;
                            if (c1 == 1)
                                tmp += (abs(hg->hash_adj[hg->size*x+y] - hgq->hash_adj[1]) + 1)*(hscrs[hgq->size*x+0] + hscrs[hgq->size*y+1]);
                            if (c2 == 1)
                                tmp += (abs(hg->hash_adj[hg->size*x+z] - hgq->hash_adj[2]) + 1)*(hscrs[hgq->size*x+0] + hscrs[hgq->size*z+2]);
                            if (c3 == 1)
                                tmp += (abs(hg->hash_adj[hg->size*y+z] - hgq->hash_adj[5]) + 1)*(hscrs[hgq->size*y+1] + hscrs[hgq->size*z+2]);

                            if (tmp < tmpm)
                            {
                                tmpm = tmp;
                                mx = x;
                                my = y;
                                mz = z;
                            }
                        }
                    }
                    if (tmpm < threshold)
                    {
                        list_cnt++;
                        tmp1->hg_node_map = (uint32_t *)malloc(3*sizeof(uint32_t));
                        tmp1->hgq_node_map = (uint32_t *)malloc(3*sizeof(uint32_t));
                        tmp1->score = (double)tmpm/((double)(c1 + c2 + c3));
                        tmp1->hg_node_map[0] = mx;
                        tmp1->hg_node_map[1] = my;
                        tmp1->hg_node_map[2] = mz;
                        tmp1->hgq_node_map[0] = 0;
                        tmp1->hgq_node_map[1] = 1;
                        tmp1->hgq_node_map[2] = 2;
                        tmp1->next = (struct motif_list *)malloc(sizeof(struct motif_list));
                        tmp2 = tmp1;
                        tmp1 = tmp1->next;
                    }
                }
            }
        }
if(list_cnt == 0){
  free(head);
  head = NULL;
}else{
        free(tmp2->next);
        tmp2->next = NULL;
}
    }
    else if (hgq->size == 4)
    {
        c1 = 0;
        c2 = 0;
        c3 = 0;
        c4 = 0;
        c5 = 0;
        c6 = 0;
        if (hgq->hash_adj[1] > 0)c1 = 1;
        if (hgq->hash_adj[2] > 0)c2 = 1;
        if (hgq->hash_adj[3] > 0)c3 = 1;
        if (hgq->hash_adj[6] > 0)c4 = 1;
        if (hgq->hash_adj[7] > 0)c5 = 1;
        if (hgq->hash_adj[11] > 0)c6 = 1;
        for (x=0; x < hg->size; x++)
        {
            for (y=0; y < hg->size; y++)
            {
                if (x != y)
                {
                    for (z=0; z < hg->size; z++)
                    {
                        if (x != z && y != z)
                        {
                            tmpm = threshold;
                            for (z1=0; z1 < hg->size; z1++)
                            {
                                if (x != z1 && y != z1 && z != z1)
                                {
                                    tmp = 0;
                                    if (c1 == 1)
                                        tmp += (abs(hg->hash_adj[hg->size*x+y] - hgq->hash_adj[1]) + 1)*(hscrs[hgq->size*x+0] + hscrs[hgq->size*y+1]);
                                    if (c2 == 1)
                                        tmp += (abs(hg->hash_adj[hg->size*x+z] - hgq->hash_adj[2]) + 1)*(hscrs[hgq->size*x+0] + hscrs[hgq->size*z+2]);
                                    if (c3 == 1)
                                        tmp += (abs(hg->hash_adj[hg->size*x+z1] - hgq->hash_adj[3]) + 1)*(hscrs[hgq->size*x+0] + hscrs[hgq->size*z1+3]);
                                    if (c4 == 1)
                                        tmp += (abs(hg->hash_adj[hg->size*y+z] - hgq->hash_adj[6]) + 1)*(hscrs[hgq->size*y+1] + hscrs[hgq->size*z+2]);
                                    if (c5 == 1)
                                        tmp += (abs(hg->hash_adj[hg->size*y+z1] - hgq->hash_adj[7]) + 1)*(hscrs[hgq->size*y+1] + hscrs[hgq->size*z1+3]);
                                    if (c6 == 1)
                                        tmp += (abs(hg->hash_adj[hg->size*z+z1] - hgq->hash_adj[11]) + 1)*(hscrs[hgq->size*z+2] + hscrs[hgq->size*z1+3]);

                                    if (tmp < tmpm)
                                    {
                                        tmpm = tmp;
                                        mx = x;
                                        my = y;
                                        mz = z;
                                        mz1 = z1;
                                    }
                                }
                            }
                            if (tmpm < threshold)
                            {
                                list_cnt++;
                                tmp1->hg_node_map = (uint32_t *)malloc(3*sizeof(uint32_t));
                                tmp1->hgq_node_map = (uint32_t *)malloc(3*sizeof(uint32_t));
                                tmp1->score = (double)tmpm/((double)(c1 + c2 + c3 + c4 + c5 + c6));
                                tmp1->hg_node_map[0] = mx;
                                tmp1->hg_node_map[1] = my;
                                tmp1->hg_node_map[2] = mz;
                                tmp1->hg_node_map[3] = mz1;
                                tmp1->hgq_node_map[0] = 0;
                                tmp1->hgq_node_map[1] = 1;
                                tmp1->hgq_node_map[2] = 2;
                                tmp1->hgq_node_map[3] = 3;
                                tmp1->next = (struct motif_list *)malloc(sizeof(struct motif_list));
                                tmp2 = tmp1;
                                tmp1 = tmp1->next;
                            }
                        }
                    }
                }
            }
        }
if(list_cnt == 0){
  free(head);
  head = NULL;
}else{
        free(tmp2->next);
        tmp2->next = NULL;
}
   }

/*

  printf("motif_list.list_cnt = %d\n\n", list_cnt);
tmp1 = head;
while(tmp1 != NULL){
  printf("%d tmp1->score = %f\n",i,tmp1->score);
  for(j = 0; j <  hgq->size; j++){
    printf("%d\t",tmp1->hg_node_map[j]);
  }
  printf("\n");
  for(j = 0; j <  hgq->size; j++){
    printf("%d\t",tmp1->hgq_node_map[j]);
  }
  i++;
  tmp1 = tmp1->next;
  printf("\n\n");
  fflush(stdout);
}

*/



  
if(list_cnt > 1){
    tmpm = head->score;
    tmp3 = head;
    tmp4 = head;
    tmp1 = head->next;
    tmp2 = head;
    while(tmp1 != NULL){
      if(tmp1->score < tmpm){
        tmpm = tmp1->score;
        tmp3 = tmp2;
        tmp4 = tmp1;
      }
      tmp2 = tmp1;
      tmp1 = tmp1->next;
    }
    if(tmp3 != tmp4){
      tmp1 = head;
      head = tmp4;
      tmp3->next = tmp4->next;
      head->next = tmp1;
    }
    tmp5 = head->next;
    tmp6 = head;
    for(x=1; x < list_cnt; x++){
      tmpm = tmp5->score;
      tmp3 = tmp5;
      tmp4 = tmp5;
      tmp1 = tmp5->next;
      tmp2 = tmp5;
      while(tmp1 != NULL){
        if(tmp1->score < tmpm){
          tmpm = tmp1->score;
          tmp3 = tmp2;
          tmp4 = tmp1;
        }
        tmp2 = tmp1;
        tmp1 = tmp1->next;
      }
      if(tmp3 != tmp4){
        tmp1 = tmp5;
        tmp5 = tmp4;
        tmp3->next = tmp4->next;
        tmp5->next = tmp1;
        tmp6->next = tmp5;
      }
      tmp6 = tmp5;
      tmp5 = tmp5->next;

    }
}


/*

  printf("motif_list.list_cnt = %d\n\n", list_cnt);
tmp1 = head;
while(tmp1 != NULL){
  printf("%d tmp1->score = %f\n",i,tmp1->score);
  for(j = 0; j <  hgq->size; j++){
    printf("%d\t",tmp1->hg_node_map[j]);
  }
  printf("\n");
  for(j = 0; j <  hgq->size; j++){
    printf("%d\t",tmp1->hgq_node_map[j]);
  }
  i++;
  tmp1 = tmp1->next;
  printf("\n\n");
  fflush(stdout);
}

*/




    ml->list_cnt = list_cnt;
    ml->next = head;
    free(hscrs);
}

void
calculateZSCORE(int subgraphCounter, uint32_t rand, int motifs, double *mean_inrand, double *var_inrand, double *total_inreal, char *path) 
{
  FILE * cm;
  int i;

  double * Score = (double *)malloc(motifs*(sizeof(double)));

  for (i = 0; i < motifs; i++) {
    mean_inrand[i] = mean_inrand[i]/rand;
    var_inrand[i] = sqrt((var_inrand[i]-(rand*(mean_inrand[i]*mean_inrand[i])))/rand);

//printf(%d m%f v%f\n",i

    if(var_inrand[i] != 0)
      Score[i] = (total_inreal[i] - mean_inrand[i])/var_inrand[i];
    else
      Score[i] = -1;
  }  

    char file[256];
    sprintf(file, "%s/ZScore.txt", path);
    printf("Writing ZScores to %s ...\n", file);
        cm = fopen(file, "w+");
    if(!cm) {
        printf("Can't open %s\n", path);
        sprintf(file, "result/ZScore.txt");
        printf("Writing to %s instead ...\n", file);
        cm = fopen(file, "w+");
        if (!cm) {
            printf("Error again ... Sorry!\n");
            exit(-1);
        }
    }

        fprintf(cm, "TOTAL NUMBER OF CLASSES:: %d\n\n", motifs);
        fprintf(cm, "ID\t\tNUM IN REAL \t\t MEAN IN RANDOM \t VAR IN RANDOM \t\t ZSCORE\n");
        for (i = 0; i < motifs; i++) {
                if (var_inrand[i] != 0)
                        fprintf(cm, "%d\t\t %f%% \t\t %f%% \t\t %f \t\t %f\n", i, total_inreal[i]/((double)subgraphCounter)*100,mean_inrand[i]/((double)subgraphCounter)*100, var_inrand[i], Score[i]);
                else
                        fprintf(cm, "**%d\t\t %f%% \t\t %f%% \t\t %f \t\t %f \n", i, total_inreal[i]/((double)subgraphCounter)*100, mean_inrand[i]/((double)subgraphCounter)*100, Score[i], (double)((total_inreal[i] - mean_inrand[i])));

        }

        fclose (cm);

  free(Score);

}


